
!+-----------------------------------------------------------------+
!PROGRAM  : 
!TYPE     : Subroutine
!PURPOSE  :   
!+-----------------------------------------------------------------+
FUNCTION random_normal() RESULT(fn_val)
  ! Adapted from the following Fortran 77 code
  !      ALGORITHM 712, COLLECTED ALGORITHMS FROM ACM.
  !      THIS WORK PUBLISHED IN TRANSACTIONS ON MATHEMATICAL SOFTWARE,
  !      VOL. 18, NO. 4, DECEMBER, 1992, PP. 434-435.
  !  The function random_normal() returns a normally distributed pseudo-random
  !  number with zero mean and unit variance.
  !  The algorithm uses the ratio of uniforms method of A.J. Kinderman
  !  and J.F. Monahan augmented with quadratic bounding curves.
  real :: fn_val
  real :: s = 0.449871, &
       t = -0.386595,   &
       a = 0.19600,     &
       b = 0.25472,     &
       r1 = 0.27597,    &
       r2 = 0.27846, u, v, x, y, q
  !Generate P = (u,v) uniform in rectangle enclosing acceptance region
  DO
     CALL RANDOM_NUMBER(u)
     CALL RANDOM_NUMBER(v)
     v = 1.7156 * (v - half)
     !     Evaluate the quadratic form
     x = u - s
     y = ABS(v) - t
     q = x**2 + y*(a*y - b*x)
     !     Accept P if inside inner ellipse
     IF (q < r1) EXIT
     !     Reject P if outside outer ellipse
     IF (q > r2) CYCLE
     !     Reject P if outside acceptance region
     IF (v**2 < -4.0*LOG(u)*u**2) EXIT
  END DO
  !     Return ratio of P's coordinates as the normal deviate
  fn_val = v/u
  RETURN
END FUNCTION random_normal
!*******************************************************************
!*******************************************************************
!*******************************************************************



!+-----------------------------------------------------------------+
!PROGRAM  : 
!TYPE     : Subroutine
!PURPOSE  :   
!+-----------------------------------------------------------------+
FUNCTION random_gamma(s, first) RESULT(fn_val)
  ! Adapted from Fortran 77 code from the book:
  !     Dagpunar, J. 'Principles of random variate generation'
  !     Clarendon Press, Oxford, 1988.   ISBN 0-19-852202-9
  !     FUNCTION GENERATES A RANDOM GAMMA VARIATE.
  !     CALLS EITHER random_gamma1 (S > 1.0)
  !     OR random_exponential (S = 1.0)
  !     OR random_gamma2 (S < 1.0).
  !     S = SHAPE PARAMETER OF DISTRIBUTION (0 < REAL).
  REAL, INTENT(IN)    :: s
  LOGICAL, INTENT(IN) :: first
  REAL                :: fn_val
  IF (s <= zero) THEN
     WRITE(*, *) 'SHAPE PARAMETER VALUE MUST BE POSITIVE'
     STOP
  END IF
  IF (s > one) THEN
     fn_val = random_gamma1(s, first)
  ELSE IF (s < one) THEN
     fn_val = random_gamma2(s, first)
  ELSE
     fn_val = random_exponential()
  END IF
  RETURN
END FUNCTION random_gamma
!*******************************************************************
!*******************************************************************
!*******************************************************************




!+-----------------------------------------------------------------+
!PROGRAM  : 
!TYPE     : Subroutine
!PURPOSE  :   
!+-----------------------------------------------------------------+

FUNCTION random_gamma1(s, first) RESULT(fn_val)
  ! Uses the algorithm in
  ! Marsaglia, G. and Tsang, W.W. (2000) `A simple method for generating
  ! gamma variables', Trans. om Math. Software (TOMS), vol.26(3), pp.363-372.
  ! Generates a random gamma deviate for shape parameter s >= 1.
  REAL, INTENT(IN)    :: s
  LOGICAL, INTENT(IN) :: first
  REAL                :: fn_val
  ! Local variables
  REAL, SAVE  :: c, d
  REAL        :: u, v, x
  IF (first) THEN
     d = s - one/3.
     c = one/SQRT(9.0*d)
  END IF
  ! Start of main loop
  DO
     ! Generate v = (1+cx)^3 where x is random normal; repeat if v <= 0.
     DO
        x = random_normal()
        v = (one + c*x)**3
        IF (v > zero) EXIT
     END DO
     ! Generate uniform variable U
     CALL RANDOM_NUMBER(u)
     IF (u < one - 0.0331*x**4) THEN
        fn_val = d*v
        EXIT
     ELSE IF (LOG(u) < half*x**2 + d*(one - v + LOG(v))) THEN
        fn_val = d*v
        EXIT
     END IF
  END DO
  RETURN
END FUNCTION random_gamma1
!*******************************************************************
!*******************************************************************
!*******************************************************************



!+-----------------------------------------------------------------+
!PROGRAM  : 
!TYPE     : Subroutine
!PURPOSE  :   
!+-----------------------------------------------------------------+
FUNCTION random_gamma2(s, first) RESULT(fn_val)
  ! Adapted from Fortran 77 code from the book:
  !     Dagpunar, J. 'Principles of random variate generation'
  !     Clarendon Press, Oxford, 1988.   ISBN 0-19-852202-9
  ! FUNCTION GENERATES A RANDOM VARIATE IN [0,INFINITY) FROM
  ! A GAMMA DISTRIBUTION WITH DENSITY PROPORTIONAL TO
  ! GAMMA2**(S-1) * EXP(-GAMMA2),
  ! USING A SWITCHING METHOD.
  !    S = SHAPE PARAMETER OF DISTRIBUTION
  !          (REAL < 1.0)
  REAL, INTENT(IN)    :: s
  LOGICAL, INTENT(IN) :: first
  REAL                :: fn_val
  !     Local variables
  REAL       :: r, x, w
  REAL, SAVE :: a, p, c, uf, vr, d
  IF (s <= zero .OR. s >= one) THEN
     WRITE(*, *) 'SHAPE PARAMETER VALUE OUTSIDE PERMITTED RANGE'
     STOP
  END IF
  IF (first) THEN                        ! Initialization, if necessary
     a = one - s
     p = a/(a + s*EXP(-a))
     IF (s < vsmall) THEN
        WRITE(*, *) 'SHAPE PARAMETER VALUE TOO SMALL'
        STOP
     END IF
     c = one/s
     uf = p*(vsmall/a)**s
     vr = one - vsmall
     d = a*LOG(a)
  END IF
  DO
     CALL RANDOM_NUMBER(r)
     IF (r >= vr) THEN
        CYCLE
     ELSE IF (r > p) THEN
        x = a - LOG((one - r)/(one - p))
        w = a*LOG(x)-d
     ELSE IF (r > uf) THEN
        x = a*(r/p)**c
        w = x
     ELSE
        fn_val = zero
        RETURN
     END IF
     CALL RANDOM_NUMBER(r)
     IF (one-r <= w .AND. r > zero) THEN
        IF (r*(w + one) >= one) CYCLE
        IF (-LOG(r) <= w) CYCLE
     END IF
     EXIT
  END DO
  fn_val = x
  RETURN
END FUNCTION random_gamma2
!*******************************************************************
!*******************************************************************
!*******************************************************************




!+-----------------------------------------------------------------+
!PROGRAM  : 
!TYPE     : Subroutine
!PURPOSE  :   
!+-----------------------------------------------------------------+
FUNCTION random_chisq(ndf, first) RESULT(fn_val)
  !     Generates a random variate from the chi-squared distribution with
  !     ndf degrees of freedom
  INTEGER, INTENT(IN) :: ndf
  LOGICAL, INTENT(IN) :: first
  REAL                :: fn_val
  fn_val = two * random_gamma(half*ndf, first)
  RETURN
END FUNCTION random_chisq

!*******************************************************************
!*******************************************************************
!*******************************************************************



!+-----------------------------------------------------------------+
!PROGRAM  : 
!TYPE     : Subroutine
!PURPOSE  :   
!+-----------------------------------------------------------------+
FUNCTION random_exponential() RESULT(fn_val)
  ! Adapted from Fortran 77 code from the book:
  !     Dagpunar, J. 'Principles of random variate generation'
  !     Clarendon Press, Oxford, 1988.   ISBN 0-19-852202-9
  ! FUNCTION GENERATES A RANDOM VARIATE IN [0,INFINITY) FROM
  ! A NEGATIVE EXPONENTIAL DlSTRIBUTION WlTH DENSITY PROPORTIONAL
  ! TO EXP(-random_exponential), USING INVERSION.
  REAL  :: fn_val
  !     Local variable
  REAL  :: r
  DO
     CALL RANDOM_NUMBER(r)
     IF (r > zero) EXIT
  END DO
  fn_val = -LOG(r)
  RETURN
END FUNCTION random_exponential
!*******************************************************************
!*******************************************************************
!*******************************************************************




!+-----------------------------------------------------------------+
!PROGRAM  : 
!TYPE     : Subroutine
!PURPOSE  :   
!+-----------------------------------------------------------------+
FUNCTION random_Weibull(a) RESULT(fn_val)
  !     Generates a random variate from the Weibull distribution with
  !     probability density:
  !                      a
  !               a-1  -x
  !     f(x) = a.x    e
  REAL, INTENT(IN) :: a
  REAL             :: fn_val
  !     For speed, there is no checking that a is not zero or very small.
  fn_val = random_exponential() ** (one/a)
  RETURN
END FUNCTION random_Weibull
!*******************************************************************
!*******************************************************************
!*******************************************************************




!+-----------------------------------------------------------------+
!PROGRAM  : 
!TYPE     : Subroutine
!PURPOSE  :   
!+-----------------------------------------------------------------+
FUNCTION random_beta(aa, bb, first) RESULT(fn_val)
  ! Adapted from Fortran 77 code from the book:
  !     Dagpunar, J. 'Principles of random variate generation'
  !     Clarendon Press, Oxford, 1988.   ISBN 0-19-852202-9
  ! FUNCTION GENERATES A RANDOM VARIATE IN [0,1]
  ! FROM A BETA DISTRIBUTION WITH DENSITY
  ! PROPORTIONAL TO BETA**(AA-1) * (1-BETA)**(BB-1).
  ! USING CHENG'S LOG LOGISTIC METHOD.
  !     AA = SHAPE PARAMETER FROM DISTRIBUTION (0 < REAL)
  !     BB = SHAPE PARAMETER FROM DISTRIBUTION (0 < REAL)
  REAL, INTENT(IN)    :: aa, bb
  LOGICAL, INTENT(IN) :: first
  REAL                :: fn_val
  !     Local variables
  REAL, PARAMETER  :: aln4 = 1.3862944
  REAL             :: a, b, g, r, s, x, y, z
  REAL, SAVE       :: d, f, h, t, c
  LOGICAL, SAVE    :: swap
  IF (aa <= zero .OR. bb <= zero) THEN
     WRITE(*, *) 'IMPERMISSIBLE SHAPE PARAMETER VALUE(S)'
     STOP
  END IF
  IF (first) THEN                        ! Initialization, if necessary
     a = aa
     b = bb
     swap = b > a
     IF (swap) THEN
        g = b
        b = a
        a = g
     END IF
     d = a/b
     f = a+b
     IF (b > one) THEN
        h = SQRT((two*a*b - f)/(f - two))
        t = one
     ELSE
        h = b
        t = one/(one + (a/(vlarge*b))**b)
     END IF
     c = a+h
  END IF
  DO
     CALL RANDOM_NUMBER(r)
     CALL RANDOM_NUMBER(x)
     s = r*r*x
     IF (r < vsmall .OR. s <= zero) CYCLE
     IF (r < t) THEN
        x = LOG(r/(one - r))/h
        y = d*EXP(x)
        z = c*x + f*LOG((one + d)/(one + y)) - aln4
        IF (s - one > z) THEN
           IF (s - s*z > one) CYCLE
           IF (LOG(s) > z) CYCLE
        END IF
        fn_val = y/(one + y)
     ELSE
        IF (4.0*s > (one + one/d)**f) CYCLE
        fn_val = one
     END IF
     EXIT
  END DO
  IF (swap) fn_val = one - fn_val
  RETURN
END FUNCTION random_beta
!*******************************************************************
!*******************************************************************
!*******************************************************************






!+-----------------------------------------------------------------+
!PROGRAM  : 
!TYPE     : Subroutine
!PURPOSE  :   
!+-----------------------------------------------------------------+
SUBROUTINE random_mvnorm(n, h, d, f, first, x, ier)
  ! Adapted from Fortran 77 code from the book:
  !     Dagpunar, J. 'Principles of random variate generation'
  !     Clarendon Press, Oxford, 1988.   ISBN 0-19-852202-9
  ! N.B. An extra argument, ier, has been added to Dagpunar's routine
  !     SUBROUTINE GENERATES AN N VARIATE RANDOM NORMAL
  !     VECTOR USING A CHOLESKY DECOMPOSITION.
  ! ARGUMENTS:
  !        N = NUMBER OF VARIATES IN VECTOR
  !           (INPUT,INTEGER >= 1)
  !     H(J) = J'TH ELEMENT OF VECTOR OF MEANS
  !           (INPUT,REAL)
  !     X(J) = J'TH ELEMENT OF DELIVERED VECTOR
  !           (OUTPUT,REAL)
  !
  !    D(J*(J-1)/2+I) = (I,J)'TH ELEMENT OF VARIANCE MATRIX (J> = I)
  !            (INPUT,REAL)
  !    F((J-1)*(2*N-J)/2+I) = (I,J)'TH ELEMENT OF LOWER TRIANGULAR
  !           DECOMPOSITION OF VARIANCE MATRIX (J <= I)
  !            (OUTPUT,REAL)
  !    FIRST = .TRUE. IF THIS IS THE FIRST CALL OF THE ROUTINE
  !    OR IF THE DISTRIBUTION HAS CHANGED SINCE THE LAST CALL OF THE ROUTINE.
  !    OTHERWISE SET TO .FALSE.
  !            (INPUT,LOGICAL)
  !    ier = 1 if the input covariance matrix is not +ve definite
  !        = 0 otherwise
  INTEGER, INTENT(IN)   :: n
  REAL, INTENT(IN)      :: h(:), d(:)   ! d(n*(n+1)/2)
  REAL, INTENT(IN OUT)  :: f(:)         ! f(n*(n+1)/2)
  REAL, INTENT(OUT)     :: x(:)
  LOGICAL, INTENT(IN)   :: first
  INTEGER, INTENT(OUT)  :: ier
  !     Local variables
  INTEGER       :: j, i, m
  REAL          :: y, v
  INTEGER, SAVE :: n2
  IF (n < 1) THEN
     WRITE(*, *) 'SIZE OF VECTOR IS NON POSITIVE'
     STOP
  END IF
  ier = 0
  IF (first) THEN                        ! Initialization, if necessary
     n2 = 2*n
     IF (d(1) < zero) THEN
        ier = 1
        RETURN
     END IF
     f(1) = SQRT(d(1))
     y = one/f(1)
     DO j = 2,n
        f(j) = d(1+j*(j-1)/2) * y
     END DO
     DO i = 2,n
        v = d(i*(i-1)/2+i)
        DO m = 1,i-1
           v = v - f((m-1)*(n2-m)/2+i)**2
        END DO
        IF (v < zero) THEN
           ier = 1
           RETURN
        END IF
        v = SQRT(v)
        y = one/v
        f((i-1)*(n2-i)/2+i) = v
        DO j = i+1,n
           v = d(j*(j-1)/2+i)
           DO m = 1,i-1
              v = v - f((m-1)*(n2-m)/2+i)*f((m-1)*(n2-m)/2 + j)
           END DO ! m = 1,i-1
           f((i-1)*(n2-i)/2 + j) = v*y
        END DO ! j = i+1,n
     END DO ! i = 2,n
  END IF
  x(1:n) = h(1:n)
  DO j = 1,n
     y = random_normal()
     DO i = j,n
        x(i) = x(i) + f((j-1)*(n2-j)/2 + i) * y
     END DO ! i = j,n
  END DO ! j = 1,n
  RETURN
END SUBROUTINE random_mvnorm
!*******************************************************************
!*******************************************************************
!*******************************************************************



!+-----------------------------------------------------------------+
!PROGRAM  : 
!TYPE     : Subroutine
!PURPOSE  :   
!+-----------------------------------------------------------------+
FUNCTION random_inv_gauss(h, b, first) RESULT(fn_val)
  ! Adapted from Fortran 77 code from the book:
  !     Dagpunar, J. 'Principles of random variate generation'
  !     Clarendon Press, Oxford, 1988.   ISBN 0-19-852202-9
  ! FUNCTION GENERATES A RANDOM VARIATE IN [0,INFINITY] FROM
  ! A REPARAMETERISED GENERALISED INVERSE GAUSSIAN (GIG) DISTRIBUTION
  ! WITH DENSITY PROPORTIONAL TO  GIG**(H-1) * EXP(-0.5*B*(GIG+1/GIG))
  ! USING A RATIO METHOD.
  !     H = PARAMETER OF DISTRIBUTION (0 <= REAL)
  !     B = PARAMETER OF DISTRIBUTION (0 < REAL)
  REAL, INTENT(IN)    :: h, b
  LOGICAL, INTENT(IN) :: first
  REAL                :: fn_val
  !     Local variables
  REAL            :: ym, xm, r, w, r1, r2, x
  REAL, SAVE      :: a, c, d, e
  REAL, PARAMETER :: quart = 0.25
  IF (h < zero .OR. b <= zero) THEN
     WRITE(*, *) 'IMPERMISSIBLE DISTRIBUTION PARAMETER VALUES'
     STOP
  END IF
  IF (first) THEN                        ! Initialization, if necessary
     IF (h > quart*b*SQRT(vlarge)) THEN
        WRITE(*, *) 'THE RATIO H:B IS TOO SMALL'
        STOP
     END IF
     e = b*b
     d = h + one
     ym = (-d + SQRT(d*d + e))/b
     IF (ym < vsmall) THEN
        WRITE(*, *) 'THE VALUE OF B IS TOO SMALL'
        STOP
     END IF
     d = h - one
     xm = (d + SQRT(d*d + e))/b
     d = half*d
     e = -quart*b
     r = xm + one/xm
     w = xm*ym
     a = w**(-half*h) * SQRT(xm/ym) * EXP(-e*(r - ym - one/ym))
     IF (a < vsmall) THEN
        WRITE(*, *) 'THE VALUE OF H IS TOO LARGE'
        STOP
     END IF
     c = -d*LOG(xm) - e*r
  END IF
  DO
     CALL RANDOM_NUMBER(r1)
     IF (r1 <= zero) CYCLE
     CALL RANDOM_NUMBER(r2)
     x = a*r2/r1
     IF (x <= zero) CYCLE
     IF (LOG(r1) < d*LOG(x) + e*(x + one/x) + c) EXIT
  END DO
  fn_val = x
  RETURN
END FUNCTION random_inv_gauss
!*******************************************************************
!*******************************************************************
!*******************************************************************



!+-----------------------------------------------------------------+
!PROGRAM  : 
!TYPE     : Subroutine
!PURPOSE  :   
!+-----------------------------------------------------------------+
FUNCTION random_Poisson(mu, first) RESULT(ival)
  !**********************************************************************
  !     Translated to Fortran 90 by Alan Miller from:
  !                           RANLIB
  !
  !     Library of Fortran Routines for Random Number Generation
  !
  !                    Compiled and Written by:
  !
  !                         Barry W. Brown
  !                          James Lovato
  !
  !             Department of Biomathematics, Box 237
  !             The University of Texas, M.D. Anderson Cancer Center
  !             1515 Holcombe Boulevard
  !             Houston, TX      77030
  !
  ! This work was supported by grant CA-16672 from the National Cancer Institute.
  !                    GENerate POIsson random deviate
  !                            Function
  ! Generates a single random deviate from a Poisson distribution with mean mu.
  !                            Arguments
  !     mu --> The mean of the Poisson distribution from which
  !            a random deviate is to be generated.
  !                              REAL mu
  !                              Method
  !     For details see:
  !               Ahrens, J.H. and Dieter, U.
  !               Computer Generation of Poisson Deviates
  !               From Modified Normal Distributions.
  !               ACM Trans. Math. Software, 8, 2
  !               (June 1982),163-179
  !     TABLES: COEFFICIENTS A0-A7 FOR STEP F. FACTORIALS FACT
  !     COEFFICIENTS A(K) - FOR PX = FK*V*V*SUM(A(K)*V**K)-DEL
  !     SEPARATION OF CASES A AND B
  !     .. Scalar Arguments ..
  REAL, INTENT(IN)    :: mu
  LOGICAL, INTENT(IN) :: first
  INTEGER             :: ival
  !     ..
  !     .. Local Scalars ..
  REAL          :: b1, b2, c, c0, c1, c2, c3, del, difmuk, e, fk, fx, fy, g,  &
       omega, px, py, t, u, v, x, xx
  REAL, SAVE    :: s, d, p, q, p0
  INTEGER       :: j, k, kflag
  LOGICAL, SAVE :: full_init
  INTEGER, SAVE :: l, m
  !     ..
  !     .. Local Arrays ..
  REAL, SAVE    :: pp(35)
  !     ..
  !     .. Data statements ..
  REAL, PARAMETER :: a0 = -.5, a1 = .3333333, a2 = -.2500068, a3 = .2000118,  &
       a4 = -.1661269, a5 = .1421878, a6 = -.1384794,   &
       a7 = .1250060
  REAL, PARAMETER :: fact(10) = (/ 1., 1., 2., 6., 24., 120., 720., 5040.,  &
       40320., 362880. /)
  !     ..
  !     .. Executable Statements ..
  IF (mu > 10.0) THEN
     !     C A S E  A. (RECALCULATION OF S, D, L IF MU HAS CHANGED)
     IF (first) THEN
        s = SQRT(mu)
        d = 6.0*mu*mu
        !             THE POISSON PROBABILITIES PK EXCEED THE DISCRETE NORMAL
        !             PROBABILITIES FK WHENEVER K >= M(MU). L=IFIX(MU-1.1484)
        !             IS AN UPPER BOUND TO M(MU) FOR ALL MU >= 10 .
        l = mu - 1.1484
        full_init = .false.
     END IF
     !     STEP N. NORMAL SAMPLE - random_normal() FOR STANDARD NORMAL DEVIATE
     g = mu + s*random_normal()
     IF (g > 0.0) THEN
        ival = g
        !     STEP I. IMMEDIATE ACCEPTANCE IF ival IS LARGE ENOUGH
        IF (ival>=l) RETURN
        !     STEP S. SQUEEZE ACCEPTANCE - SAMPLE U
        fk = ival
        difmuk = mu - fk
        CALL RANDOM_NUMBER(u)
        IF (d*u >= difmuk*difmuk*difmuk) RETURN
     END IF
     !     STEP P. PREPARATIONS FOR STEPS Q AND H.
     !             (RECALCULATIONS OF PARAMETERS IF NECESSARY)
     !             .3989423=(2*PI)**(-.5)  .416667E-1=1./24.  .1428571=1./7.
     !             THE QUANTITIES B1, B2, C3, C2, C1, C0 ARE FOR THE HERMITE
     !             APPROXIMATIONS TO THE DISCRETE NORMAL PROBABILITIES FK.
     !             C=.1069/MU GUARANTEES MAJORIZATION BY THE 'HAT'-FUNCTION.
     IF (.NOT. full_init) THEN
        omega = .3989423/s
        b1 = .4166667E-1/mu
        b2 = .3*b1*b1
        c3 = .1428571*b1*b2
        c2 = b2 - 15.*c3
        c1 = b1 - 6.*b2 + 45.*c3
        c0 = 1. - b1 + 3.*b2 - 15.*c3
        c = .1069/mu
        full_init = .true.
     END IF
     IF (g < 0.0) GO TO 50
     !             'SUBROUTINE' F IS CALLED (KFLAG=0 FOR CORRECT RETURN)
     kflag = 0
     GO TO 70
     !     STEP Q. QUOTIENT ACCEPTANCE (RARE CASE)
40   IF (fy-u*fy <= py*EXP(px-fx)) RETURN
     !     STEP E. EXPONENTIAL SAMPLE - random_exponential() FOR STANDARD EXPONENTIAL
     !             DEVIATE E AND SAMPLE T FROM THE LAPLACE 'HAT'
     !             (IF T <= -.6744 THEN PK < FK FOR ALL MU >= 10.)
50   e = random_exponential()
     CALL RANDOM_NUMBER(u)
     u = u + u - one
     t = 1.8 + SIGN(e, u)
     IF (t <= (-.6744)) GO TO 50
     ival = mu + s*t
     fk = ival
     difmuk = mu - fk
     !             'SUBROUTINE' F IS CALLED (KFLAG=1 FOR CORRECT RETURN)
     kflag = 1
     GO TO 70
     !     STEP H. HAT ACCEPTANCE (E IS REPEATED ON REJECTION)\
60   IF (c*ABS(u) > py*EXP(px+e) - fy*EXP(fx+e)) GO TO 50
     RETURN
     !     STEP F. 'SUBROUTINE' F. CALCULATION OF PX, PY, FX, FY.
     !             CASE ival < 10 USES FACTORIALS FROM TABLE FACT
70   IF (ival>=10) GO TO 80
     px = -mu
     py = mu**ival/fact(ival+1)
     GO TO 110
     !             CASE ival >= 10 USES POLYNOMIAL APPROXIMATION
     !             A0-A7 FOR ACCURACY WHEN ADVISABLE
     !             .8333333E-1=1./12.  .3989423=(2*PI)**(-.5)
80   del = .8333333E-1/fk
     del = del - 4.8*del*del*del
     v = difmuk/fk
     IF (ABS(v)>0.25) THEN
        px = fk*LOG(one + v) - difmuk - del
     ELSE
        px = fk*v*v* (((((((a7*v+a6)*v+a5)*v+a4)*v+a3)*v+a2)*v+a1)*v+a0) - del
     END IF
     py = .3989423/SQRT(fk)
110  x = (half - difmuk)/s
     xx = x*x
     fx = -half*xx
     fy = omega* (((c3*xx + c2)*xx + c1)*xx + c0)
     IF (kflag <= 0) GO TO 40
     GO TO 60
     !---------------------------------------------------------------------------
     !     C A S E  B.    mu < 10
     !     START NEW TABLE AND CALCULATE P0 IF NECESSARY
  ELSE
     IF (first) THEN
        m = MAX(1, INT(mu))
        l = 0
        p = EXP(-mu)
        q = p
        p0 = p
     END IF
     !     STEP U. UNIFORM SAMPLE FOR INVERSION METHOD
     DO
        CALL RANDOM_NUMBER(u)
        ival = 0
        IF (u <= p0) RETURN
        !     STEP T. TABLE COMPARISON UNTIL THE END PP(L) OF THE
        !             PP-TABLE OF CUMULATIVE POISSON PROBABILITIES
        !             (0.458=PP(9) FOR MU=10)
        IF (l == 0) GO TO 150
        j = 1
        IF (u > 0.458) j = MIN(l, m)
        DO k = j, l
           IF (u <= pp(k)) GO TO 180
        END DO
        IF (l == 35) CYCLE
        !     STEP C. CREATION OF NEW POISSON PROBABILITIES P
        !             AND THEIR CUMULATIVES Q=PP(K)
150     l = l + 1
        DO k = l, 35
           p = p*mu / k
           q = q + p
           pp(k) = q
           IF (u <= q) GO TO 170
        END DO
        l = 35
     END DO
170  l = k
180  ival = k
     RETURN
  END IF
  RETURN
END FUNCTION random_Poisson
!*******************************************************************
!*******************************************************************
!*******************************************************************



!+-----------------------------------------------------------------+
!PROGRAM  : 
!TYPE     : Subroutine
!PURPOSE  :   
!+-----------------------------------------------------------------+  
FUNCTION random_binomial1(n, p, first) RESULT(ival)
  ! FUNCTION GENERATES A RANDOM BINOMIAL VARIATE USING C.D.Kemp's method.
  ! This algorithm is suitable when many random variates are required
  ! with the SAME parameter values for n & p.
  !    P = BERNOULLI SUCCESS PROBABILITY
  !           (0 <= REAL <= 1)
  !    N = NUMBER OF BERNOULLI TRIALS
  !           (1 <= INTEGER)
  !    FIRST = .TRUE. for the first call using the current parameter values
  !          = .FALSE. if the values of (n,p) are unchanged from last call
  ! Reference: Kemp, C.D. (1986). `A modal method for generating binomial
  !            variables', Commun. Statist. - Theor. Meth. 15(3), 805-813.
  INTEGER, INTENT(IN) :: n
  REAL, INTENT(IN)    :: p
  LOGICAL, INTENT(IN) :: first
  INTEGER             :: ival
  !     Local variables
  INTEGER         :: ru, rd
  INTEGER, SAVE   :: r0
  REAL            :: u, pd, pu
  REAL, SAVE      :: odds_ratio, p_r
  REAL, PARAMETER :: zero = 0.0, one = 1.0
  IF (first) THEN
     r0 = (n+1)*p
     p_r = bin_prob(n, p, r0)
     odds_ratio = p / (one - p)
  END IF
  CALL RANDOM_NUMBER(u)
  u = u - p_r
  IF (u < zero) THEN
     ival = r0
     RETURN
  END IF
  pu = p_r
  ru = r0
  pd = p_r
  rd = r0
  DO
     rd = rd - 1
     IF (rd >= 0) THEN
        pd = pd * (rd+1) / (odds_ratio * (n-rd))
        u = u - pd
        IF (u < zero) THEN
           ival = rd
           RETURN
        END IF
     END IF
     ru = ru + 1
     IF (ru <= n) THEN
        pu = pu * (n-ru+1) * odds_ratio / ru
        u = u - pu
        IF (u < zero) THEN
           ival = ru
           RETURN
        END IF
     END IF
  END DO
  !     This point should not be reached, but just in case:
  ival = r0
  RETURN
END FUNCTION random_binomial1
!*******************************************************************
!*******************************************************************
!*******************************************************************





!+-----------------------------------------------------------------+
!PROGRAM  : 
!TYPE     : Subroutine
!PURPOSE  :   
!+-----------------------------------------------------------------+
FUNCTION bin_prob(n, p, r) RESULT(fn_val)
  !     Calculate a binomial probability
  INTEGER, INTENT(IN) :: n, r
  REAL, INTENT(IN)    :: p
  REAL                :: fn_val
  !     Local variable
  REAL                :: one = 1.0
  fn_val = EXP( lngamma(DBLE(n+1)) - lngamma(DBLE(r+1)) - lngamma(DBLE(n-r+1)) &
       + r*LOG(p) + (n-r)*LOG(one - p) )
  RETURN
END FUNCTION bin_prob
!*******************************************************************
!*******************************************************************
!*******************************************************************




!+-----------------------------------------------------------------+
!PROGRAM  : 
!TYPE     : Subroutine
!PURPOSE  :   
!+-----------------------------------------------------------------+
FUNCTION lngamma(x) RESULT(fn_val)
  ! Logarithm to base e of the gamma function.
  !
  ! Accurate to about 1.e-14.
  ! Programmer: Alan Miller
  ! Latest revision of Fortran 77 version - 28 February 1988
  REAL (dp), INTENT(IN) :: x
  REAL (dp)             :: fn_val
  !       Local variables
  REAL (dp) :: a1 = -4.166666666554424D-02, a2 = 2.430554511376954D-03,  &
       a3 = -7.685928044064347D-04, a4 = 5.660478426014386D-04,  &
       temp, arg, product, lnrt2pi = 9.189385332046727D-1,       &
       pi = 3.141592653589793D0
  LOGICAL   :: reflect
  !       lngamma is not defined if x = 0 or a negative integer.
  IF (x > 0.d0) GO TO 10
  IF (x /= INT(x)) GO TO 10
  fn_val = 0.d0
  RETURN
  !       If x < 0, use the reflection formula:
  !               gamma(x) * gamma(1-x) = pi * cosec(pi.x)
10 reflect = (x < 0.d0)
  IF (reflect) THEN
     arg = 1.d0 - x
  ELSE
     arg = x
  END IF
  !       Increase the argument, if necessary, to make it > 10.
  product = 1.d0
20 IF (arg <= 10.d0) THEN
     product = product * arg
     arg = arg + 1.d0
     GO TO 20
  END IF
  !  Use a polynomial approximation to Stirling's formula.
  !  N.B. The real Stirling's formula is used here, not the simpler, but less
  !       accurate formula given by De Moivre in a letter to Stirling, which
  !       is the one usually quoted.
  arg = arg - 0.5D0
  temp = 1.d0/arg**2
  fn_val = lnrt2pi + arg * (LOG(arg) - 1.d0 + &
       (((a4*temp + a3)*temp + a2)*temp + a1)*temp) - LOG(product)
  IF (reflect) THEN
     temp = SIN(pi * x)
     fn_val = LOG(pi/temp) - fn_val
  END IF
  RETURN
END FUNCTION lngamma
!*******************************************************************
!*******************************************************************
!*******************************************************************




!+-----------------------------------------------------------------+
!PROGRAM  : 
!TYPE     : Subroutine
!PURPOSE  :   
!+-----------------------------------------------------------------+
FUNCTION random_binomial2(n, pp, first) RESULT(ival)
  !**********************************************************************
  !     Translated to Fortran 90 by Alan Miller from:
  !                              RANLIB
  !
  !     Library of Fortran Routines for Random Number Generation
  !
  !                      Compiled and Written by:
  !
  !                           Barry W. Brown
  !                            James Lovato
  !
  !               Department of Biomathematics, Box 237
  !               The University of Texas, M.D. Anderson Cancer Center
  !               1515 Holcombe Boulevard
  !               Houston, TX      77030
  !
  ! This work was supported by grant CA-16672 from the National Cancer Institute.
  !                    GENerate BINomial random deviate
  !                              Function
  !     Generates a single random deviate from a binomial
  !     distribution whose number of trials is N and whose
  !     probability of an event in each trial is P.
  !                              Arguments
  !     N  --> The number of trials in the binomial distribution
  !            from which a random deviate is to be generated.
  !                              INTEGER N
  !     P  --> The probability of an event in each trial of the
  !            binomial distribution from which a random deviate
  !            is to be generated.
  !                              REAL P
  !     FIRST --> Set FIRST = .TRUE. for the first call to perform initialization
  !               the set FIRST = .FALSE. for further calls using the same pair
  !               of parameter values (N, P).
  !                              LOGICAL FIRST
  !     random_binomial2 <-- A random deviate yielding the number of events
  !                from N independent trials, each of which has
  !                a probability of event P.
  !                              INTEGER random_binomial
  !                              Method
  !     This is algorithm BTPE from:
  !         Kachitvichyanukul, V. and Schmeiser, B. W.
  !         Binomial Random Variate Generation.
  !         Communications of the ACM, 31, 2 (February, 1988) 216.
  !**********************************************************************
  !*****DETERMINE APPROPRIATE ALGORITHM AND WHETHER SETUP IS NECESSARY
  !     ..
  !     .. Scalar Arguments ..
  REAL, INTENT(IN)    :: pp
  INTEGER, INTENT(IN) :: n
  LOGICAL, INTENT(IN) :: first
  INTEGER             :: ival
  !     ..
  !     .. Local Scalars ..
  REAL            :: alv, amaxp, f, f1, f2, u, v, w, w2, x, x1, x2, ynorm, z, z2
  REAL, PARAMETER :: zero = 0.0, half = 0.5, one = 1.0
  INTEGER         :: i, ix, ix1, k, mp
  INTEGER, SAVE   :: m
  REAL, SAVE      :: p, q, xnp, ffm, fm, xnpq, p1, xm, xl, xr, c, al, xll,  &
       xlr, p2, p3, p4, qn, r, g
  !     ..
  !     .. Executable Statements ..
  !*****SETUP, PERFORM ONLY WHEN PARAMETERS CHANGE
  IF (first) THEN
     p = MIN(pp, one-pp)
     q = one - p
     xnp = n * p
  END IF
  IF (xnp > 30.) THEN
     IF (first) THEN
        ffm = xnp + p
        m = ffm
        fm = m
        xnpq = xnp * q
        p1 = INT(2.195*SQRT(xnpq) - 4.6*q) + half
        xm = fm + half
        xl = xm - p1
        xr = xm + p1
        c = 0.134 + 20.5 / (15.3 + fm)
        al = (ffm-xl) / (ffm - xl*p)
        xll = al * (one + half*al)
        al = (xr - ffm) / (xr*q)
        xlr = al * (one + half*al)
        p2 = p1 * (one + c + c)
        p3 = p2 + c / xll
        p4 = p3 + c / xlr
     END IF
     !*****GENERATE VARIATE, Binomial mean at least 30.
20   CALL RANDOM_NUMBER(u)
     u = u * p4
     CALL RANDOM_NUMBER(v)
     !     TRIANGULAR REGION
     IF (u <= p1) THEN
        ix = xm - p1 * v + u
        GO TO 110
     END IF
     !     PARALLELOGRAM REGION
     IF (u <= p2) THEN
        x = xl + (u-p1) / c
        v = v * c + one - ABS(xm-x) / p1
        IF (v > one .OR. v <= zero) GO TO 20
        ix = x
     ELSE
        !     LEFT TAIL
        IF (u <= p3) THEN
           ix = xl + LOG(v) / xll
           IF (ix < 0) GO TO 20
           v = v * (u-p2) * xll
        ELSE
           !     RIGHT TAIL
           ix = xr - LOG(v) / xlr
           IF (ix > n) GO TO 20
           v = v * (u-p3) * xlr
        END IF
     END IF
     !*****DETERMINE APPROPRIATE WAY TO PERFORM ACCEPT/REJECT TEST
     k = ABS(ix-m)
     IF (k <= 20 .OR. k >= xnpq/2-1) THEN
        !     EXPLICIT EVALUATION
        f = one
        r = p / q
        g = (n+1) * r
        IF (m < ix) THEN
           mp = m + 1
           DO i = mp, ix
              f = f * (g/i-r)
           END DO
        ELSE IF (m > ix) THEN
           ix1 = ix + 1
           DO i = ix1, m
              f = f / (g/i-r)
           END DO
        END IF
        IF (v > f) THEN
           GO TO 20
        ELSE
           GO TO 110
        END IF
     END IF
     !     SQUEEZING USING UPPER AND LOWER BOUNDS ON LOG(F(X))
     amaxp = (k/xnpq) * ((k*(k/3. + .625) + .1666666666666)/xnpq + half)
     ynorm = -k * k / (2.*xnpq)
     alv = LOG(v)
     IF (alv<ynorm - amaxp) GO TO 110
     IF (alv>ynorm + amaxp) GO TO 20
     !     STIRLING'S (actually de Moivre's) FORMULA TO MACHINE ACCURACY FOR
     !     THE FINAL ACCEPTANCE/REJECTION TEST
     x1 = ix + 1
     f1 = fm + one
     z = n + 1 - fm
     w = n - ix + one
     z2 = z * z
     x2 = x1 * x1
     f2 = f1 * f1
     w2 = w * w
     IF (alv - (xm*LOG(f1/x1) + (n-m+half)*LOG(z/w) + (ix-m)*LOG(w*p/(x1*q)) +    &
          (13860.-(462.-(132.-(99.-140./f2)/f2)/f2)/f2)/f1/166320. +               &
          (13860.-(462.-(132.-(99.-140./z2)/z2)/z2)/z2)/z/166320. +                &
          (13860.-(462.-(132.-(99.-140./x2)/x2)/x2)/x2)/x1/166320. +               &
          (13860.-(462.-(132.-(99.-140./w2)/w2)/w2)/w2)/w/166320.) > zero) THEN
        GO TO 20
     ELSE
        GO TO 110
     END IF
  ELSE
     !     INVERSE CDF LOGIC FOR MEAN LESS THAN 30
     IF (first) THEN
        qn = q ** n
        r = p / q
        g = r * (n+1)
     END IF
90   ix = 0
     f = qn
     CALL RANDOM_NUMBER(u)
100  IF (u >= f) THEN
        IF (ix > 110) GO TO 90
        u = u - f
        ix = ix + 1
        f = f * (g/ix - r)
        GO TO 100
     END IF
  END IF
110 IF (pp > half) ix = n - ix
  ival = ix
  RETURN
END FUNCTION random_binomial2
!*******************************************************************
!*******************************************************************
!*******************************************************************




!+-----------------------------------------------------------------+
!PROGRAM  : 
!TYPE     : Subroutine
!PURPOSE  :   
!+-----------------------------------------------------------------+
FUNCTION random_neg_binomial(sk, p) RESULT(ival)
  ! Adapted from Fortran 77 code from the book:
  !     Dagpunar, J. 'Principles of random variate generation'
  !     Clarendon Press, Oxford, 1988.   ISBN 0-19-852202-9
  ! FUNCTION GENERATES A RANDOM NEGATIVE BINOMIAL VARIATE USING UNSTORED
  ! INVERSION AND/OR THE REPRODUCTIVE PROPERTY.
  !    SK = NUMBER OF FAILURES REQUIRED (Dagpunar's words!)
  !       = the `power' parameter of the negative binomial
  !           (0 < REAL)
  !    P = BERNOULLI SUCCESS PROBABILITY
  !           (0 < REAL < 1)
  ! THE PARAMETER H IS SET SO THAT UNSTORED INVERSION ONLY IS USED WHEN P <= H,
  ! OTHERWISE A COMBINATION OF UNSTORED INVERSION AND
  ! THE REPRODUCTIVE PROPERTY IS USED.
  REAL, INTENT(IN)   :: sk, p
  INTEGER            :: ival
  !     Local variables
  ! THE PARAMETER ULN = -LOG(MACHINE'S SMALLEST REAL NUMBER).
  REAL, PARAMETER    :: h = 0.7
  REAL               :: q, x, st, uln, v, r, s, y, g
  INTEGER            :: k, i, n
  IF (sk <= zero .OR. p <= zero .OR. p >= one) THEN
     WRITE(*, *) 'IMPERMISSIBLE DISTRIBUTION PARAMETER VALUES'
     STOP
  END IF
  q = one - p
  x = zero
  st = sk
  IF (p > h) THEN
     v = one/LOG(p)
     k = st
     DO i = 1,k
        DO
           CALL RANDOM_NUMBER(r)
           IF (r > zero) EXIT
        END DO
        n = v*LOG(r)
        x = x + n
     END DO
     st = st - k
  END IF
  s = zero
  uln = -LOG(vsmall)
  IF (st > -uln/LOG(q)) THEN
     WRITE(*, *) ' P IS TOO LARGE FOR THIS VALUE OF SK'
     STOP
  END IF
  y = q**st
  g = st
  CALL RANDOM_NUMBER(r)
  DO
     IF (y > r) EXIT
     r = r - y
     s = s + one
     y = y*p*g/s
     g = g + one
  END DO
  ival = x + s + half
  RETURN
END FUNCTION random_neg_binomial
!*******************************************************************
!*******************************************************************
!*******************************************************************



!+-----------------------------------------------------------------+
!PROGRAM  : 
!TYPE     : Subroutine
!PURPOSE  :   
!+-----------------------------------------------------------------+
FUNCTION random_von_Mises(k, first) RESULT(fn_val)
  !     Algorithm VMD from:
  !     Dagpunar, J.S. (1990) `Sampling from the von Mises distribution via a
  !     comparison of random numbers', J. of Appl. Statist., 17, 165-168.
  !     Fortran 90 code by Alan Miller
  !     CSIRO Division of Mathematical & Information Sciences
  !     Arguments:
  !     k (real)        parameter of the von Mises distribution.
  !     first (logical) set to .TRUE. the first time that the function
  !                     is called, or the first time with a new value
  !                     for k.   When first = .TRUE., the function sets
  !                     up starting values and may be very much slower.
  REAL, INTENT(IN)     :: k
  LOGICAL, INTENT(IN)  :: first
  REAL                 :: fn_val
  !     Local variables
  INTEGER          :: j, n
  INTEGER, SAVE    :: nk
  REAL, PARAMETER  :: pi = 3.14159265
  REAL, SAVE       :: p(20), theta(0:20)
  REAL             :: sump, r, th, lambda, rlast
  REAL (dp)        :: dk
  IF (first) THEN                        ! Initialization, if necessary
     IF (k < zero) THEN
        WRITE(*, *) '** Error: argument k for random_von_Mises = ', k
        RETURN
     END IF
     nk = k + k + one
     IF (nk > 20) THEN
        WRITE(*, *) '** Error: argument k for random_von_Mises = ', k
        RETURN
     END IF
     dk = k
     theta(0) = zero
     IF (k > half) THEN
        !     Set up array p of probabilities.
        sump = zero
        DO j = 1, nk
           IF (j < nk) THEN
              theta(j) = ACOS(one - j/k)
           ELSE
              theta(nk) = pi
           END IF
           !     Numerical integration of e^[k.cos(x)] from theta(j-1) to theta(j)
           CALL integral(theta(j-1), theta(j), p(j), dk)
           sump = sump + p(j)
        END DO
        p(1:nk) = p(1:nk) / sump
     ELSE
        p(1) = one
        theta(1) = pi
     END IF                         ! if k > 0.5
  END IF                           ! if first
  CALL RANDOM_NUMBER(r)
  DO j = 1, nk
     r = r - p(j)
     IF (r < zero) EXIT
  END DO
  r = -r/p(j)
  DO
     th = theta(j-1) + r*(theta(j) - theta(j-1))
     lambda = k - j + one - k*COS(th)
     n = 1
     rlast = lambda
     DO
        CALL RANDOM_NUMBER(r)
        IF (r > rlast) EXIT
        n = n + 1
        rlast = r
     END DO
     IF (n .NE. 2*(n/2)) EXIT         ! is n even?
     CALL RANDOM_NUMBER(r)
  END DO
  fn_val = SIGN(th, (r - rlast)/(one - rlast) - half)
  RETURN
END FUNCTION random_von_Mises
!*******************************************************************
!*******************************************************************
!*******************************************************************



!+-----------------------------------------------------------------+
!PROGRAM  : 
!TYPE     : Subroutine
!PURPOSE  :   
!+-----------------------------------------------------------------+
SUBROUTINE integral(a, b, result, dk)
  !     Gaussian integration of exp(k.cosx) from a to b.
  REAL (dp), INTENT(IN) :: dk
  REAL, INTENT(IN)      :: a, b
  REAL, INTENT(OUT)     :: result
  !     Local variables
  REAL (dp)  :: xmid, range, x1, x2,                                    &
       x(3) = (/0.238619186083197_dp, 0.661209386466265_dp, 0.932469514203152_dp/), &
       w(3) = (/0.467913934572691_dp, 0.360761573048139_dp, 0.171324492379170_dp/)
  INTEGER    :: i
  xmid = (a + b)/2._dp
  range = (b - a)/2._dp
  result = 0._dp
  DO i = 1, 3
     x1 = xmid + x(i)*range
     x2 = xmid - x(i)*range
     result = result + w(i)*(EXP(dk*COS(x1)) + EXP(dk*COS(x2)))
  END DO
  result = result * range
  RETURN
END SUBROUTINE integral
!*******************************************************************
!*******************************************************************
!*******************************************************************



!+-----------------------------------------------------------------+
!PROGRAM  : 
!TYPE     : Subroutine
!PURPOSE  :   
!+-----------------------------------------------------------------+
FUNCTION random_Cauchy() RESULT(fn_val)
  !     Generate a random deviate from the standard Cauchy distribution
  REAL     :: fn_val
  !     Local variables
  REAL     :: v(2)
  DO
     CALL RANDOM_NUMBER(v)
     v = two*(v - half)
     IF (ABS(v(2)) < vsmall) CYCLE               ! Test for zero
     IF (v(1)**2 + v(2)**2 < one) EXIT
  END DO
  fn_val = v(1) / v(2)
  RETURN
END FUNCTION random_Cauchy
!*******************************************************************
!*******************************************************************
!*******************************************************************
